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1. Introduction 

It is possible to estimate the mass and radius of a transiting planet using a combination 
of photometry, Doppler data, and stellar modeling, just as has long been done for eclipsing 
binary stars (Vogcl 1890). Some aspects of the planetary orbit and atmosphere can also be 
measured through high-precision transit photometry and spectroscopy (see, e.g., Charbon- 
neau et al. 2007, Seager 2008, Winn 2008). These opportunities are more likely to occur 
for short-period planets, because the probability for a randomly-oriented orbit to be viewed 
close enough to edge-on for transits declines as P~^/^. This explains why all but one of the 
known transiting planets have periods smaller than 10 days. 

The exception is HD 17156b, for which P — 21.2 days. This planet was discovered 
in a Doppler survey by Fischer et al. (2007). The probability for transits to occur was 

larger than one might have guessed based only on the period, because the planet is near 
pericenter at the time of inferior conjunction (see, e.g.. Burke 2008, Barnes 2008, or Kane 
& von Braun 2008), and indeed transits were discovered by Barbieri et al. (2007). The 
relatively long period, along with the large orbital eccentricity of 0.67, presents an interesting 
opportunity to study the planetary atmospheric response to strongly time-variable heating 
by the parent star (Langton & Laughlin 2007). However, the long period also presents a 
challenge to observers. Prom a given site, there are only 2-3 good opportunities each year to 
observe a complete transit. Irwin et al. (2008), Narita et al. (2008), and Gillon et al. (2008) 
observed the photometric transits of Oct. 1, Nov. 12, and Dec. 3, 2007, respectively. They 
were able to improve upon the precision of the system parameters given originally by Barbieri 
et al. (2007), but left further scope for improvement through higher- precision photometry. 

In this paper we present photometry of the transit of Dec. 25, 2007 based on observations 
with 4 different telescopes. We also present 10 new measurements of the Doppler shift of the 
host star, gathered outside of transits. We analyze all of these data to refine the estimates of 
the system parameters, using similar techniques to those we have applied previously as part 
of the Transit Light Curve (TLC) project (Holman et al. 2006, Winn et al. 2007) and that 
were recently applied to a sample of 23 planets by Torres et al. (2008). The observations 
and data reduction are described in § 2. The data analysis is described in § 3. In § 4, 
the results of the data analysis are used together with other observed stellar properties and 
stellar evolutionary models to determine the properties of the star and planet. All of the 
results are summarized in § 5, and potential future studies are discussed. 
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2. Observations and Data Reduction 

2.1. Out-of-transit radial velocities 

Fischer et al. (2007) reported 33 measurements of the radial velocity (RV) of HD 17156 
over a time range from January 2006 to February 2007. Of these, 9 velocities were obtained 
with the Subaru 8m telescope and HDS spectrograph, with a precision of approximately 
5 m s~^. The other 24 velocities were based on observations with the Keck 1 10m telescope 
and HIRES spectrograph (Vogt et al. 1994), with a precision of 1-2 m s~^. 

To these, we add 10 new Keck/HIRES velocities with 1-2 m s~^ precision that were 
obtained between August 2007 and March 2008. The new data are based on observations with 
the same telescope, instrument, and setup as the previous Keck observations. In particular we 
employed an iodine gas absorption cell to calibrate the instrumental profile and wavelength 
scale. To maintain a consistent signal-to-noise ratio of about 200 per resolution element, we 
employed the HIRES exposure meter, which uses a pickoff mirror to direct a small fraction 
of the starhght to a photomultipher tube and monitors the exposure level as a function of 
time (Kibrick et al. 2006). 

All 34 of the Keck/HIRES velocities were re- measured based on an improved reduction of 
the raw CCD images and a refined version of the algorithm of Butler et al. (1996), including 
the use of a new stellar template. (The stellar template, an important ingredient in the 
deconvolution algorithm of Butler et al. 1996, is an observation of the target star obtained 
without the iodine cell and with a higher signal-to-noise ratio and higher resolution than 
the rest of the spectra.) The resulting velocities are given in Tabled! For convenience, this 
table also includes the Subaru/HDS velocities reported previously. Fig. [1] shows the radial 
velocities as a function of orbital phase, using the ephemeris derived in § 3. 

2.2. Transit photometry 

We observed the transit of UT 25 December 2007 (JD 2454459) using telescopes at two 
different observatories in Arizona: the Fred L. Whipple Observatory (FLWO) and Fairborn 
Observatory. The Moon was full. The weather was generally clear over both observatories, 
although there were some light clouds and transparency variations. 

At FLWO, we used the 1.2m telescope and Keplercam, a 4096^ CCD with a 23' x 23' 
field of view (Szentgyorgyi et al. 2005). The images were binned 2x2, giving a scale of 0'.'68 
per binned pixel. We obtained 10 s exposures though a 2;-band filter for 6.5 hr bracketing 
the predicted midtransit time. The telescope was defocused to avoid significant nonlinearity 
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Fig. 1. — Radial velocity (RV) variation of HD 17156. Top. — The measured RVs as a 
function of orbital phase, expressed in days since midtransit. Bottom. — Differences between 
the observed and calculated RVs, using the model described in § 3.3. The root-mean-squared 
(RMS) residual is 3.8 m s'^. 
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Table 1. Relative Radial Velocities (RV) of HD 17156 



Observatory Code'' Heliocentric Julian Date RV [m s ^] Measurement Uncertainty [m s ^] 
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^(1) Keck/HIRES. (2) Subaru/HDS. 

Note. — The time stamps represent the Heliocentric Julian Date at the time of midexposure. 
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and saturation, giving a typical full width at half-maximum (FWHM) of stellar images of 
4" (6 pixels). The time between exposures was 11 s due to readout and reset operations. 
Autoguiding failures led to pointing drifts of 10 pixels in right ascension and 40 pixels in 
declination over the course of the night. During the observations the target rose from airmass 
1.34 to 1.30 (reaching the meridian at HJD 2454459.67) and then set to airmass 1.64. 

We used standard IRAF procedures for overscan correction, trimming, bias subtraction, 
flat-field division, and aperture photometry of HD 17156 and 32 other stars in the field of 
view (all necessarily fainter than HD 17156). We created a reference signal by combining 
the normalized light curves of the different comparison stars, and divided the flux history of 
HD 17156 by this reference signal. We experimented with different choices for the aperture 
size, combinations of the comparison stars, and weighting schemes for combining the nor- 
malized light curves, aiming to minimize the standard deviation of the out-of-transit portion 
of the light curve. Best results were obtained when the aperture diameter was 15 pixels and 
the comparison light curves were weighted by the square root of the mean flux. 

At Fairborn Observatory, we used 3 independent telescopes: the T8, TIO, and Til 0.8m 
automated photometric telescopes (APTs). All of the APTs are equipped with two temperature- 
stabilized EMI 9124QB photomultiplier tubes for measuring photon count rates simultane- 
ously through Stromgren b and y filters. Each telescope nodded back and forth between 
HD 17156 {V = 8.17, B -V = 0.64) and the comparison star HD 15784 (V = 6.64, 
B — V = 0.64), which was found to be constant in brightness by Fischer et al. (2007). The 
T8, TIO, and Til APTs obtained 170, 190, and 200 observations of both stars, respectively, 
over a period 7.2 hours. From these measurements, we computed a total of 560 target-minus- 
comparison differential magnitudes in each b and y photometric band. The time series was 
trimmed to 523 points to eliminate bad data taken at the highest airmass. To increase the 
signal-to-noise ratio in the resulting light curves, we averaged the b and y passbands together, 
resulting in a synthetic (6 -|- y)/2 passband. 

The FLWO and APT light curves were corrected for differing airmass extinction be- 
tween the target star and comparison stars, by fitting the observed magnitudes to a linear 
function of airmass (see Eqn. [3]). The parameters of the linear function were determined si- 
multaneously with the other model parameters, as described in § 3. The extinction-corrected 
data are given in Table [21 and plotted in Fig. [21 along with the best-fitting model. For the 
APT data, the standard deviation of the out-of-transit (OOT) data and of the residuals are 
0.00214 and 0.00213, respectively. The corresponding figures for the FLWO data are 0.00233 
and 0.00206, respectively. For both light curves, the expected photometric precision due 
only to Poisson noise and scintillation noise (using the empirical formulas of Reiger 1963, 
Young 1967 and Dravins 1998) is about 50% smaller than the OOT standard deviation, with 
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Fig. 2.— Relative photometry of HD 17156 on UT 25 Dec 2007, along with the best-fitting 
model. Top. — Data from the T8, TIO, and Til 0.8m automated photometric telescopes at 
Fairborn Observatory, in the "(fe + |/)/2" band (the Stromgren b and y data were averaged). 
The median time interval between data points is 43 s. The standard deviation of the out- 
of-transit (OOT) data and of the residuals are 0.00214 and 0.00213, respectively. Bottom. — 
Data from the Fred L. Whipple 1.2m telescope and Keplercam, in the z band. The median 
time interval between data points is 24 s. The standard deviation of the OOT data and of 
the residuals are 0.00233 and 0.00206, respectively. 
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the dominant contribution arising from scintillation noise. It is possible that the scintillation 
noise was higher than the empirical and highly approximate formula predicts. It is also likely 
that there are additional noise sources such as transparency variations beyond the simple 
dependence on airmass. 



At the heart of our analysis was a simultaneous fit of a parametric model to the new 
light curves and the available radial- velocity data. In § 3.1 we describe the model and 
fitting procedures. In § 3.2 we provide more detail about the photometric noise and how it 
was modeled. In § 3.3 we discuss the determination of midtransit times and a new transit 
ephemeris. In § 4 we explain how the data analysis was integrated with stellar evolutionary 
models to determine the final system parameters. 



The model is based on a two-body Keplerian orbit, with the loss of light during tran- 
sits given by the formulae of Mandel & Agol (2002). The orbit is parameterized by the 
semi-amplitude (K), period (P), midtransit time (Tcjll, inclination (i), eccentricity (e), and 
argument of pericenter (uj). The coordinate system is chosen such that i < 90° and the 
longitude of nodes is zero. 

Additional parameters relevant to the photometric data are the planet-to-star radius ra- 
tio (Rp/R^,), the semimajor axis in units of the stellar radius [a/R^), and the limb-darkening 
coefficients {ui and U2) for each light curve. We assumed a quadratic limb- darkening law, 



where n is the cosine of the angle between the line of sight and the normal to the stellar 
surface, and is the specific intensity. In many studies of transiting planets, the limb- 
darkening coefficients are often held fixed at values deemed appropriate for the host star. 



^Since the orbit of HD 17156b is highly eccentric, the definition of Tc requires some care. Here we define 
it as the time when the projected planet-star separation is smallest, which is also the time of minimum light 
for a limb-darkened star. This is to be distinguished from other possible definitions, such as the halfway 
point between first and last contact, or the moment when the true anomaly / of the planetary orbit is equal 
to 7r/2 — uj. 



3. Analysis of radial-velocity and photometric data 



3.1. Joint radial velocity and photometric analysis 
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based on stellar- atmosphere models. However, when the data are sufficiently precise it is 
preferable to fit for the coefficients, as recently emphasized by Southworth (2008). This 
is because the model atmospheres might be wrong, and at least are uncertain to some 
degree. Neglecting this uncertainty leads to underestimated errors in all parameters that are 
covariant with the limb-darkening coefficients, namely Rp/R^, R^/a, and i. For this study 
we allowed ui and U2 to vary freely subject only to the conditions ui + U2 < I (nonnegative 
intensity at the limb), Ui + U2 > (fainter at the limb than the center), and t^i > 
(maximum intensity at the center). It proved advantageous to perform the fit using the 
linear combinations 

7 7 

Vl^Ui + -U2, V2 = --Ml + U2, (2) 

because vi and V2 have nearly uncorrelated errors (for further discussion, see Pal 2008). 

The model also has 6 nuisance parameters. For each of the two RV data sets there is an 
additive constant velocity (7^ for Keck/HIRES and 75 for Subaru/HDS). For each of the 
two light curves there are the parameters mo and k of the differential extinction correction, 

mcor = mobs + mo + kz, (3) 

where mobs is the observed magnitude, mcor is the corrected magnitude, and z is the airmass. 

All together there were 18 parameters: K, P, Tc, i, e, uj, Rp/R^, a/Ri,^ an additive 
velocity for each of two RV data sets, two limb-darkening coefficients Ui and U2 for each of 
two bandpasses, and two parameters itlq and k for each of two bandpasses to describe the 
correction for differential airmass extinction. The fitting statistic was 

= + + Xl, (4) 
with the terms defined as follows. The first term is based on the fit to the photometric data: 



i=l 



^^^^ '"/i(obs) - /i(calc) 



(5) 



where /j(obs) is the ith measured flux (photometric data point), /j(calc) is the calculated 
flux given a particular choice of model parameters, and af^i is the uncertainty in the ith 
measured flux. The choice of cx/^j is nontrivial and is discussed in § 3.2. Here we simply state 
our choice to take cr/^j to be a constant for each light curve given by cr/^j = [3a^es) where a^^s 
is the standard deviation of the flux residuals, and /3 > 1 is a factor intended to account 
for time-correlated errors. For the APT light curve, a^es — 0.00213 and (3 — 1.03, giving 
an effective uncertainty per point of 0.0022. For the FLWO fight curve, Ures = 0.00207 and 
(3 — 1.28, giving an effective uncertainty per point of 0.0026. 
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The second term in Eqn. (jlj) is based on the fit to the RV data: 

2 v-^ r^(obs) — V^(calc)l ^ 




(6) 



where l^(obs) and l^(calc) are the ith observed and calculated RVs, and av,i is the corre- 
sponding uncertainty. For (Ty^i, we used the quadrature sum of the measurement uncertainty 
given in Table [1] and a constant cxyo = 3.4 m s~^ intended to account for stellar "jitter," 
excess noise that is usually attributed to motions of the stellar photosphere. This choice 
resulted in a reduced Xv unity when fitting the RV data only, and it is consistent with 
observations of other stars of similar type (Wright 2005). 

The third term in Eqn. (jl]) is an a priori constraint enforcing the transit ephemeris that 
is derived in § 3.3. Specifically, 



The central values of P and Tc, and their uncertainties, are derived in § 3.3 based on ob- 
servations of 5 different transits spanning 106 days. They are more precise than could be 
derived internally from only the RV data, the APT data, and the FLWO data that are fitted 
here; hence during this step we treated P and Tc as externally measured quantities. Since 
the transit ephemeris is based in part on the midtransit times of the APT and FLWO light 
curves, we used an iterative procedure: first, a previously published ephemeris was used in 
Eqn. ([7j); second, the midtransit times based on the APT and FLWO light curves were de- 
termined as in § 3.3; and third, the ephemeris was re-derived and used in the next iteration 
of the fit to the RV and photometric data. Further iterations were performed but made no 
appreciable difference in the results. 

We used a Markov Chain Monte Carlo (MCMC) algorithm to estimate the best-fitting 
model parameters and their uncertainties (see, e.g.. Appendix A of Tegmark et al. 2004). 
This algorithm creates a sequence of points (a "chain") in parameter space by iterating 
a jump function, which in our case was the addition of a Gaussian random deviate to a 
randomly-selected single parameter. After this operation, if the new point has a lower 
than the previous point, the "jump" is executed: the new point is added to the chain. If not, 
then the jump is executed with a probability proportional to exp(— Ax^/2). If the jump is 
not executed, the current point is repeated in the chain. The sizes of the random deviates 
are set to values for which ~40% of jumps are executed. After creating multiple chains to 
check for mutual convergence, and trimming off the initial segments to eliminate artifacts of 
the initial condition, the density of the chain's points in parameter space is taken to be the 
joint a posteriori probability distribution of the parameter values. Probability distributions 




P (days) -21.21688 
0.00044 



+ 



Tc (HJD) - 2, 454, 459.69987 
0.00045 



(7) 
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for individual parameters are created by marginalizing over all other parameters. For each 
parameter we report the mode of the distribution and the 68.3% confidence limits, defined 
by the 15.85% percentile and the 84.15% percentile in the cumulative distribution. 

The results are given in Table [3l with the designation A in the last column. The 
entries designated B are those that are drawn from other works and are repeated here for 
convenience. The entries designated C are based on a synthesis of our modeling results and 
theoretical models of stellar evolution, as discussed in § 4. The results for the limb-darkening 
parameters are given separately, in Table |H 



Deriving reliable uncertainties in the system parameters requires a realistic treatment 
of the photometric noise. In this work, as in others, we have based our analysis on a 
statistic that implicitly treats the noise as independent and identically-distributed Gaussian 
random variables added to each datum. This is because simple and fast statistical methods 
are applicable in that case, and because there is no clear alternative. Unfortunately, precise 
photometry is often plagued with time-correlated ( "red" ) noise, which is probably responsible 
for the frequent disagreements between reported values of photometric parameters in the 
literature. It is advisable to attempt to justify the common assumption of independent 
Gaussian noise, or at least to gain some understanding of the limitations of that assumption]^. 

To investigate the noise in our light curves, we performed some tests on the residuals. 
The top two panels in Fig. [3] show the residuals for the two independent light curves. The 
second row shows histograms of the residuals, which are approximately Gaussian. The third 
row shows the autocorrelation of the residuals as a function of the lag. 



where N is the number of data points, is the ith residual flux (measured flux minus 
calculated flux) and f is the mean residual flux. The autocorrelations are ^0.1. The fourth 
row shows the Allan deviation o"^(/) of the residuals, where a a is defined by (Allan 1964; 



^In this vein we quote an email from G. Kovacs: "Identifying the 'red noise component' (whatever it 
means) in a single realization of a rather limited time series (such as photometric foUowup data) is a mission 
impossible." 



3.2. Photometric noise analysis 





Thompson, Moran, & Swensen 2001) 

2 

(9) 

where N is the number of data points, denotes the ith residual, and I is the lag. The 
Allan deviation is commonly used in the time metrology literature to assess 1// noise. The 
solid lines show the dependence aA oc that is expected of white Gaussian noise. The 
data fall close to these lines. The bottom row of panels shows the standard deviation of the 
time-binned residuals, as a function of the number of data points per bin. Specifically, we 
averaged the residuals into m bins of n points and calculated the standard deviation (T„ of 
the binned residuals. In the absence of red noise, one would expect 

(10) 

77,^/^ \m — 1/ 

but in reality an is larger than this expression by a factor (3. In general (3 depends on m but 
the dependence is weak. For bin sizes ranging from 10-30 min (the most important timescale 
of the transit light curve) the mean value of /3 is 1.03 for the APT light curve and 1.28 for 
the FLWO light curve. Apparently the APT residuals average down as one would expect of 
white noise, while the FLWO residuals are correlated to some degree. 

As stated above, our x^-based method is strictly appropriate only for uncorrelated 
errors. One might imagine modifying the definition of Xf employ the full covariance 
matrix of the errors rather than assuming independent errors. However, given that there 
is little structure in the autocorrelation functions to provide guidance on how to model the 
covariance matrix, and that the correlations seem small, we account for correlated noise 
in a simple and approximate fashion: we assign a photometric error bar of cr/^j = /JcTres to 
each data point, where a^es is the standard deviation of the residuals. The choice /3 — 1 
is essentially equivalent the common procedure of assuming equal errors for all data points 
and scaling the error bars such that x^/A^dof = 1- We used /3 = 1.03 for the APT data and 
P = 1.28 for the FLWO data. 



1 



2{N+l-l) 



N-21 
n=0 



1 ^ 



m=0 



3.3. Midtransit times 

Measurements of midtransit times are important for refining the transit ephemeris and 
thereby enabling accurate predictions for future observations, and also for searching for 
satellites and additional planets via the method of Holman & Murray (2005) and Agol et 
al. (2005). We determined the midtransit time from each light curve using 3 different methods 
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Fig. 3. — Photometric noise analysis. The left column applies to the APT {b + y)/2 data, 
and the right column applies to the FLWO 1.2m z data. First row. — Residuals from the 
best-fitting model (observed— calculated). Second row. — Histogram of the residuals. Third 
row. — Autocorrelation of the residuals as a function of the lag. Fourth row. — Allan deviation 
of the residuals, as a function of the lag. The straight line shows the l'^^"^ dependence that is 
expected of white Gaussian noise. Fifth row. — Standard deviation of time-binned residuals, 
as a function of the number of data points per bin. The straight line shows the expectation 
for white Gaussian noise (see Eqn. [TOj) . 
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to check for consistency: the MCMC algorithm described in § 3.1, and two different bootstrap 
analyses described below. All of these methods rely on the photometric transit model and 
the computation of Xf- Besides Tc, the only other variable parameters were the slope (k) and 
offset (mo) of the differential extinction correction. These were the only other parameters 
covariant with the midtransit time. All other parameters needed to specify the transit model 
were held fixed at the best-fitting values. 

The first bootstrap method was "random draws with replacement" (Press et al. 1992, 
p. 689). It involves minimizing Xf a function of the parameters for 10^ synthetic data sets, 
each of which has the same number of data points as the real data. Each entry in a synthetic 
data set is a datum (a time stamp and relative flux) drawn randomly from the real data 
set, with repetitions allowed. Thus, a substantial fraction of the entries in each synthetic 
data set are duplicated at least once and receive greater weight in Xf sum. The idea is to 
estimate the noise properties of the data using the observed data values themselves, rather 
than choosing models for the underlying physical process and for the noise. However, an 
underlying assumption is that the errors are uncorrelated and identically distributed. The 
distribution of results for Tc is taken to be the probability density for the midtransit time. 

The second bootstrap method was "residual permutation." It is similar to the method 
just described but the synthetic data sets are created differently. The residuals of the best- 
fitting model are added back to the model light curve after performing a cyclic permutation 
of their time indices. With N data points, one may create — 1 synthetic data sets in this 
manner. Another A^ — 1 may be created by inverting the time order of the residuals and then 
performing cyclic permutations. The idea is again to use the data themselves to estimate 
the noise properties, but this time without assuming that the errors are uncorrelated. The 
correlations between residuals at different times are preserved. The distribution of results 
for Tc is taken to be the probability density for the midtransit time. A disadvantage of this 
method is that only 2(A^ — 1) realizations can be generated without further assumptions, 
and therefore the distribution is relatively noisy. 

The results from all 3 methods are given in Table [51 As before, the quoted values are 
the modes of the probability distributions and the quoted error bars range from the 15.85% 
percentile to the 84.15% percentile. The distributions are symmetric in all cases, and are 
nearly Gaussian for the MCMC and random-draws methods. The residual-permutation 
method produced distributions with broader wings than a Gaussian function. For each light 
curve, the results from all 3 methods are in agreement within 0.15a, where a is the error 
in the MCMC method. The bootstrap method gave the smallest error bars, as one might 
expect, given that the bootstrap method ignores correlated noise. (The bootstrap-derived 
error bar is smaller than the MCMC-derived error bar by approximately the "red noise" 
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factor j3.) In what follows, we adopt the MCMC results for concreteness and for consistency 
with our previous analyses. The error bars on the midtransit times derived from the FLWO 
and APT light curves are 1.1 min and 0.9 min, and the difference between the results is 
1.0 min. This level of agreement is a consistency check on the accuracy of our error bars. 

We fitted a linear ephemeris, Tc[E] = Tc[0] + EP, to all of the midtransit times at our 
disposal: namely, the APT and FLWO midtransit times presented in this paper, and the 4 
different midtransit times reported by Barbieri et al. (2007), Gillon et al. (2008), Irwin et 
al. (2008), and Narita et al. (2008). For convenience, all of the midtransit times are given in 
Table [61 The results were 

Tc[0] = 2, 454, 459.69987 ± 0.00045 [HJD], (11) 
P = 21.21688 ± 0.00044 days. (12) 

The fit gives = 1-99 with 4 degrees of freedom, suggesting that a constant period is 
consistent with the available data. A plot of the timing residuals (observed— calculated 
midtransit times) is shown in Fig. HI There are no obvious anomalies at the level of a few 
minutes. 

I I I I I I _ 



-6 -5 -4 -3 -2-10 1 

Epoch 

Fig. 4. — Transit timing residuals for HD 17156b. The calculated times, using the ephemeris 
given in Eqns. fllllll2p . have been subtracted from the observed times. 
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4. Theoretical isochrone fitting 

The RV and photometric data do not uniquely determine the masses and radii of the 
planet and the star. Some external information about the star or the planet must be intro- 
duced to break the fitting degeneracies Mp oc mI'^ and Rp (x oc mI'^ (see, e.g., Winn 
2008). We broke these degeneracies by requiring consistency between the observed properties 
of the star, the stellar mean density that can be derived from the photometric parameter 
a/Ri, (Seager & Mallen-Ornelas 2003, Sozzetti et al. 2007), and theoretical models of stellar 
evolution. The inputs were Teg = 6079 ± 80 K and [Fe/H] = +0.24 ± 0.05 (from Fischer et 
al. 2007 but with enlarged error bars, as per Torres et al. 2008), the absolute magnitud^ 
My = 3.80 ± 0.12, the stellar mean density p^ = 0.589+°;°!;^ g cm derived from the results 
for the a/ Ri, parameter, and the Yonsei-Yale (Y^) stellar evolution models by Yi et al. (2001) 
and Demarque et al. (2004). We computed isochrones for the allowed range of metallicities, 
and for stellar ages ranging from 0.1 to 14 Gyr. For each stellar property (mass, radius, and 
age), we took a weighted average of the points on each isochrone, in which the weights were 
proportional to exp(— Ax^/2) with 
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Here, the A quantities denote the deviations between the observed and calculated values 
at each point. The asymmetric error bar in p^ was taken into account by using different 
values of depending on the sign of the deviation. The weights were further multiplied by 
a factor taking into account the number density of stars along each isochrone, assuming a 
Salpeter mass function. We used the same code as Torres et al. (2008) and refer the reader 
to that paper for further details. 

Through this analysis, we found = 1.263+EJ;|]^^ Mq, R^ = LUGtof^r ^0. and a stellar 
age of 3.18^Q gg Gyr. The corresponding planetary mass and radius were obtained by merging 
the results for the stellar properties with the parameters determined in our analysis of the 
RV and photometric data. The results are Mp = 3.2121°;°^^ Mj^p and Rp = 1.023l|j:[}^g R^^p. 
These values are also given in Table [31 along with the values for some other interesting 
parameters that can be derived from the preceding results. 

As a consistency check, we computed the implied stellar surface gravity and its un- 
certainty based on our analysis, finding logg = 4.219^o;q55 where g is in cm s~^. This 



■^The quoted My is based on the transformation of Tycho-2 apparent magnitudes (H0g et al. 2000) to 
the Johnson V band, giving V = 8.172 ± 0.012, and the Hipparcos parallax of tt = 13.34 ± 0.72 mas (van 
Leeuwen 2007). We have assumed zero reddening. 
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agrees with the spectroscopic determination of surface gravity, logg = 4.29 ± 0.06 (Fischer 
et ah 2007), based on an analysis of the widths of pressure-sensitive hues in the optical 
spectrum. Since the error bars on the two results are comparable, one might be tempted 
to use the spectroscopic logs' ^ further constraint on the stellar properties. We did not 
take this approach out of concern that the spectroscopic determination is more complex and 
liable to underestimation of the error (see, e.g., Winn et al. 2008). 

5. Discussion 

5.1. Summary and comparison with previous results 

We have presented photometry of a transit of HD 17156b, and 10 additional measure- 
ments of the radial velocity of the host star. We have analyzed the data along with the 
other observed stellar properties to refine the estimates of the basic system parameters. For 
many of the parameters, our results are more precise than those reported previously, despite 
some aspects of the previous analyses that may have caused the quoted error bars to be un- 
realistically small. Specifically, none of the previous analyses allowed for any uncertainty in 
the limb-darkening law, and neither Barbieri et al. (2007) nor Irwin et al. (2008) attempted 
to quantify the effect of time-correlated noise. Another improvement in our analysis was 
the integration of the data analysis with stellar evolution models to make full use of the 
information in the light curve and arrive at a self-consistent solution. 

Our result for the stellar mass, = 1.263lQ Q4y agrees with the previous de- 
termination of 1.2 ± 0.1 Mq by Fischer et al. (2007) and improves on the precision. Our 
result for the stellar radius, 1.446]'^g;og7 Rq, agrees with the Fischer et al. (2007) estimate of 
1.470 ± 0.085 Rq. The essential difference between the two analyses is that we made use 
of the photometric determination of p^, while Fischer et al. (2007) used the spectroscopic 
determination of logs'. As for the planet, our result of 1.0231q;q55 -Rjup is very similar to 
the value 1.01 ± 0.09 -Rjup found by Irwin et al. (2008). This should be interpreted as an 
agreement between the measured transit depths, and not necessarily the other light curve 
parameters, because Irwin et al. (2008) did not fit for the a/R^, parameter. Instead they used 
an external constraint on that parameter based on the work by Fischer et al. (2007). Gillon 
et al. (2008) also found the same transit depth, but a smaller value of a/-R* and therefore 
larger values of R^, and Rp (by about la, where a is the error quoted by Gillon et al. 2008). 
Our result for Rp is approximately 3 times more precise. 

The results for the limb-darkening coefficients, given in Table IH show that the darkening 
is greater in the {b + y)/2 band than in the z band, as expected. The values of the coefficients 
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themselves are poorly constrained and the errors are highly correlated, making it difficult 
to compare to specific theoretical predictions. For the z band, interpolation of the tables of 
Claret (2004) gives a prediction of U\ = 0.17 and U2 = 0.35, which is just outside the la 
range of our results. The data favor a less limb-darkened star. The tables of Claret (2000) 
give Ml = 0.52 and U2 = 0.28 for the b band, and ui = 0.39 and U2 = 0.33 for the y band. 
Either set of coefficients is compatible with the loose bounds provided by the APT data. 



5.2. Comparison with theoretical models 

A primary goal of precise transit observations is to compare the observed planetary 
properties with theoretical models of the planet's interior structure. For example, a persistent 
theme in this field is that at least a few planets have radii that are "too large" by the 
standards of theoretical models of solar-composition giant planets, even after accounting 
for the intense stellar heating and selection effects (see Burrows et al. 2007 for a recent 
discussion). Other planets are so small that the models fit only when the composition is 
altered to be much richer in heavy elements than the Sun, a possible indication of the dense 
interior cores that are expected according to the core-accretion theory of planet formation 
(see, e.g., Sato et al. 2005). 

Bodenheimer et al. (2003) give predictions for the radii of giant planets as a function of 
the age, mass, and time-averaged equilibrium temperature of the planet, defined as 



T 



= (783 K) (1 - Ay/\ (14) 



167raa2(l + eV2)2 



where A is the Bond albedo, is the stellar luminosity, a is the Stefan-Boltzmann constant, 
a is the semimajor axis, and e is the eccentricity. In the latter equality we have evaluated 
Teq for HD 17156b using the results given in Table [31 As long as the albedo is not very 
close to unity, Bodenheimer et al. (2003) predict a planetary radius of 1.10 i?jup for a solar 
composition at 4.5 Gyr. They also considered enriching the solar composition by 40 of 
additional heavy elements (4% of the total mass) and found that the radius decreases by less 
than 1%. 

Fortney et al. (2007) have presented theoretical models parameterized by mass, age, 
and an effective orbital distance, defined as the distance from the Sun where a hypothetical 
planet on a circular orbit would receive the same time- averaged fiux as the actual planet, 

/r \i/2 

d^ = a{l + e^/2) i-^j = 0.13 AU, (15) 
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where again we have evaluated the expression as appropriate for HD 17156b. Fortney et 
al. (2007) also predicted 1.10 Rj^p for a solar composition at 4.5 Gyr, decreasing to 1.02 i?jup 
when 100 of solar-composition material (10% of the total mass) is replaced by a heavy- 
element "core." This would provide a good match to the observed radius of 1.023 Rjup- 

We conclude from these comparisons, as did Irwin et al. (2008), that HD 17156b is 
smaller than a theoretical giant planet of solar-composition, and that this may be an in- 
dication of heavy-element enrichment. Heavy-element enrichment is expected according to 
the core-accretion model of planet formation (Mizuno 1980, Pollack et al. 1996). There is 
evidence for such enrichment in Jupiter and Saturn, with Jupiter in particular consisting 
of 3-15% heavy elements (Guillot 2005). One should be wary of over-interpretation, given 
that our measurement of Rp differs from 1.1 Rj^p by only 1.1 times the measurement un- 
certainty. However, there are some known factors that would increase the theoretical radius 
and thereby enlarge the discrepancy at least slightly: the fiducial radius calculated by Bo- 
denheimer et al. (2003) and Fortney et al. (2007) refers to a higher pressure (smaller radius) 
than the transit-measured radius (Burrows et al. 2003); the models do not take into account 
tidal heating and consequent inflation due to the nonzero eccentricity (Liu et al. 2008); and 
the age of the system is estimated to be 3 Gyr, younger than the 4.5 Gyr age for which the 
models were calculated. 



5.3. The spin-orbit angle, and the probability of secondary eclipses 

One parameter that we have not improved on, but that deserves mention, is the angle 
between the stellar spin axis and the orbital axis. One might expect the axes to be well- 
aligned, given that the star and planet formed from a common disk and given the good 
alignment observed in the Solar system. Then again, by the same logic the orbit should be 
circular, and it is not. A lower limit on the angle between the spin and orbital axes can 
be derived by monitoring the apparent Doppler shift throughout a transit, exploiting the 
Rossiter- McLaughlin (RM) effect (see, e.g., Queloz et al. 2000, Winn et al. 2005). Using 
this technique, Narita et al. (2008) found that the angle between sky projections of the two 
axes is A = 62 ± 25 deg, a 2.5o" misalignment. Cochran et al. (2008) presented two different 
data sets giving A = 4.5 ± 15.6 deg and A = —32.4 ± 25.2 deg, and concluded that the 
data are consistent with good alignment. Unfortunately, our new data alone do not allow 
for significant progress on this issue, because the limiting errors arc in the precision of the 
RM data, which we have not improved in this study. Further spectroscopic observations are 
warranted. 

Another important angle is the orbital inclination with respect to the sky plane. If the 
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orbit is oriented close enough to edge-on, then the planet will be periodically eclipsed by 
the parent star. Observations of such secondary eclipses would reveal the planetary albedo 
or thermal emission from the planet, depending upon the observing bandpass. In addition, 
infrared observations could help to understand the radiative dynamics of the planetary at- 
mosphere, as emphasized by Barbieri et al. (2007). To check whether secondary eclipses 
are likely to occur, we used our MCMC results to compute the a posteriori probability 
distribution for the impact parameter at superior conjunction. 



We find the probability for secondary eclipses {bu < 1 + Rp/Ri,) to be 18%. The probability 
for complete eclipses, or occultations (6n < 1 + Rp/R*), is 15%. These odds are better than 
those found by Irwin et al. (2008), which were 6.9% and 9.2%. Gillon et al. (2008) found the 
probability of occultations to be 0.04%. Presumably this significant difference is attributable 
to our finding of a more edge-on orbit (i = 86.2tQ;8 deg, as opposed to 85.5t deg from Gillon 
et al. (2008). Our error bars for i are no smaller than those reported previously because the 
achievable error in i worsens rapidly as i approaches 90 deg (Carter et al. 2008). To help 
in planning observations, we have used our results to predict the timing of the events, as 
well as the quantity {Rp/duY (where du is the star-planet distance at superior conjunction) 
which sets the amplitude of the reflected-light signal from the planet. The results for this 
latter parameter are conditioned on the assumption that secondary eclipses do indeed occur. 
The results are given in Table [71 With a bit of luck, HD 17156b will be eclipsed and give 
observers another gift. 
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Table 2. Relative Photometry of HD 17156 



Observatory Code* Heliocentric Julian Date Relative flux 



2454459 


56340 


1 


0041 


2454459 


56360 


1 


0019 


2454459 


56490 


1 


0017 


2454459 


56500 


1 


0006 


2454459 


56630 





9954 


2454459 


56650 


1 


0005 


2454459 


56800 


1 


0016 


2454459 


56920 





9964 


2454459 


57060 





9974 


2454459 


57090 


1 


0018 


2454459 


57210 


1 


0009 


2454459 


57230 


1 


0011 



''(1) T8, TIO, and Til APT 0.8m telescopes, Fairborn Observa- 
tory, Arizona, USA. (2) Fred L. Whipple Observatory 1.2m tele- 
scope, Arizona, USA. 

Note. — The time stamps represent the Heliocentric Julian Date 
at the time of mid-cxposurc. Wc intend for this Table to appear in 
entirety in the electronic version of the journal. An excerpt is shown 
here to illustrate its format. The data are also available from the 
authors upon request. 
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Tablc 3. System Parameters of HD 17156 



Parameter Value 68.3% Conf. Limits Comment 



Transit and orbital parameters: 



Orbital period, P [d] 


21.21688 


±0.00044 


A 


Midtransit time [HJD] 


2454459.69987 


±0.00045 


A 


Planct-to-star radius ratio, Rp/Ri, 


0.0727 


±0.0016 


A 


Orbital inclination, i [dcg] 


86.2 


-0.8, +2.1 


A 


Scaled semimajor axis, a/B,,, 


22.8 


-1.7, +2.2 


A 


Transit impact parameter, bi 


0.55 


-0.29, +0.03 


A 


Transit duration [hr] 


3.177 


-0.041, +0.071 


A 


Transit ingress or egress duration [hr] 


0.25 


-0.026, +0.078 


A 


Velocity semiamplitude, K [m s^^] 


272.7 


±2.1 


A 


Orbital eccentricity, e 


0.675.3 


±0.0036 


A 


Argument of pericenter, to [deg] 


121.64 


±0.48 


A 


Planet-to-star mass ratio, M^jMi, 


0.00244 


±0.000029 


C 


Semimajor axis [AU] 


0.1623 


-0.0020, +0.0015 


C 


Stellar parameters: 








Mass, Mi, [Mq] 


1.263 


-0.047, +0.035 


C 


Radius, R« [R©] 


1.446 


-0.067, +0.099 


C 


Surface gravity, log [cgs] 


4.219 


-0.055, +0.033 





Mean density, p* [g cm~*] 


0.589 


-0.103, +0.066 


A 


Effective temperature, Tgfi [K] 


6079 


±80 


B 


Metallicity, [Fc/H] 


-i-0.24 


±0.05 


B 


Projected rotation rate, vsiaii, [km s~^] 


2.6 


±0.5 


B 


Luminosity [L©] 


2.55 


-0.32, +0.24 


C 


Absolute V magnitude 


3.78 


-0.27,+0.23 


C 


Age [Gyr] 


3.06 


-0.76, +0.64 


C 


Planetary parameters: 








Mp [Mjup] 


3.212 


-0.082,+0.069 


C 


Rp [Fi„7up] 


1.023 


-0.055,+0.070 


C 


Surface gravity, gp [m s~'^] 


76.1 


-9.0,+7.3 


A 


Mean density, pp [g cm~^] 


3.72 


±0.67 


C 



Note. — (A) Based on the joint analysis of photometric and RV data (see § 3.1-3.3). (B) From 
Fischer et al. (2007), with enlarged error bars for T^g and [Fe/H]. (0) Functions of group A and 
B parameters, supplemented as needed by an isochrone analysis (see § 4), the Tycho-2 apparent 
magnitudes (H0g et al. 2000), and Hipparcos parallax of tt = 13.34 ±0.72 mas (van Leeuwen 2007). 



Table 4. Fitted Limb-Darkening Coefficients for HD 17156 



Bandpass 


Linear Coefficient (u\) 


Quadratic Coefficient (^2) 


z 

(6 + !/)/2 


'-'•^'-'-0.06 


o.ool°;?| 

0.64l°:i? 
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Table 5. Comparison of Methods for Measuring Midtransit Times 



Analysis method Midtransit time [HJD] 

FLWO APT 



MCMC with /3 > 1 2454159.70044 ± 0.00075 2454459.69972 ± 0.00065 

Random draws witii replacement 2451459.70034 ± 0.00054 2454459.69974 ± 0.00063 
Residual permulal ion 2 15 1 159.700 18 ± 0,00090 2 15 1 15<).0<)!)81 ± 0,00095 



Table 6. Midtransit Times of HD 17156b 



Heliocentric Julian Date 


1(7 Uncertainty 


Reference 


2454353.61000 


0.02000 


Barbieri et al. (2007) 


2454438.48271 


0.00067 


Gillon ct al. (2008) 


2454374.83380 


0.00200 


Irwin et al. (2008) 


2454417.26450 


0.00210 


Narita et al. (2008) 


2454459.70044 


0.00075 


This work (FLWO) 


2454459.69972 


0.00065 


This work (APT) 



Table 7. Predicted Superior Conjunction Parameters for HD 17156b 



Parameter 


Value 


68.3% Conf. Limits 


Midpoint of superior conjunction [HJD] 


2454464.627 


-0.090, +0.100 


Probability of occultation, P(6ii < 1 + Rp/R*) 


18% 




Probability of non-grazing occultation, P{bii < 1 — Rp/R*) 


15% 




Occultation duration for an edge-on orbit [hr] 


12.9 


-0.9,+1.4 


Reflected-light figure of merit, 10^ (Rp/du)^ 


4.68 


-0.22, +0.31 



Note. — The midpoint of superior conjunction has not been corrected for the light-travel time 
across the system. The result for (Rp/du)^ is conditioned on the occurrence of occultations. 



